###R Code for Figure 1
#install.packages("gamlss")
#install.packages("XLConnect")

setwd("F:/Count Modeling Replication/zigerell/")
library(gamlss)
library(gamlss.dist)
library(XLConnect)

#Read in 'unchangedoriginal.xlsx' from working directory
zigerell<-readWorksheetFromFile('unchangedoriginaldata.xlsx',sheet=1,header=TRUE)

graphics.off()



###Run newfunc.R first

##Negative Binomial Graph
mNBI <- histDist1(zigerell$sscie_count, "NBI", main = "NBI", trace = FALSE)

tt1<-mNBI$tt1  		## access observed frequencies
tt2<-c(tt1,rep(0,1100-length(tt1)))	## augment some zeros to the end to make it have a "nice" length
names(tt2)<-c(names(tt1),as.character((length(tt1)+1):length(tt2)))		##give augmented bits names

yy1<-mNBI$yy1						## access fitted probabilities
yy2<-c(yy1,rep(0,1100-length(yy1)))	

group<-rep(1:(1100/50),each=50)			## create some groups each of length 50 (this can be anything as long as it divides 1100)

tt3<-tapply(tt2,group,sum)						## sum up observed frequencies in each group		
yy3<-tapply(yy2,group,sum)						## sum up fitted probabilities
names(tt3)<-tapply(as.numeric(names(tt2)),group,mean)		## 

r<-barplot(tt3,ylim = c(0,1),, fg = "blue", 
           col = "gray", axis.lty = 1, border = "blue", 
           col.axis = "blue4", col.main = "blue4", col.lab = "blue4", main="Fitted NB Distribution",
           ylab="Frequency Density/Fitted Probabilities", xlab="Counts")
lines(r, yy3, type = "h", col = "red", lwd = 2, col.axis = "blue", 
      col.main = "blue4", col.lab = "blue4", fg = gray(0.7))
points(r, yy3, col = "red", pch = 21, lwd = 2, col.axis = "blue")



#Poisson Inverse Gaussian Graph
mPIG <- histDist1(zigerell$sscie_count, "PIG", main = "PIG", trace = FALSE)

tt1<-mPIG$tt1      ## access observed frequencies
tt2<-c(tt1,rep(0,1100-length(tt1)))	## augment some zeros to the end to make it have a "nice" length
names(tt2)<-c(names(tt1),as.character((length(tt1)+1):length(tt2)))		##give augmented bits names

yy1<-mPIG$yy1						## access fitted probabilities
yy2<-c(yy1,rep(0,1100-length(yy1)))	

group<-rep(1:(1100/50),each=50)			## create some groups each of length 50 (this can be anything as long as it divides 1100)

tt3<-tapply(tt2,group,sum)						## sum up observed frequencies in each group		
yy3<-tapply(yy2,group,sum)						## sum up fitted probabilities
names(tt3)<-tapply(as.numeric(names(tt2)),group,mean)		## 

r<-barplot(tt3,ylim = c(0,1),, fg = "blue", 
           col = "gray", axis.lty = 1, border = "blue", 
           col.axis = "blue4", col.main = "blue4", col.lab = "blue4", main="Fitted PIG Distribution",
           ylab="Frequency Density/Fitted Probabilities", xlab="Counts")
lines(r, yy3, type = "h", col = "red", lwd = 2, col.axis = "blue", 
      col.main = "blue4", col.lab = "blue4", fg = gray(0.7))
points(r, yy3, col = "red", pch = 21, lwd = 2, col.axis = "blue")
